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Nitrogen isotope analysis of bone collagen has been used to reconstruct the breastfeeding practices of 
archaeological human populations. However, weaning ages have been estimated subjectively because 
of a lack of both information on subadult bone collagen turnover rates and appropriate analytical 
models. Here, we present a model for analyzing cross-sectional 5 15 l\l data of subadult bone collagen, 
which incorporates newly estimated bone collagen turnover rates and a framework of approximate 
Bayesian computation. Temporal changes in human subadult bone collagen turnover rates were 
pq \ estimated anew from data on tissue-level bone metabolism reported in previous studies. A model 

for reconstructing precise weaning ages was then developed and incorporating the estimated turnover 
rates. The model is presented as a new open source R package, WARN (Weaning Age Reconstruction 
with Nitrogen isotope analysis), which computes the age at the start and end of weaning, 15 N- 
enrichment through maternal to infant tissue, and <5 15 N value of collagen synthesized entirely from 



weaning foods with their posterior probabilities. A precise reconstruction of past breastfeeding and 
weaning practices over a wide range of time periods and geographic regions could make it possible to 
understand this unique feature of human life history and cultural diversity in infant feeding practices. 



Introduction 

Investigating variations in the breastfeeding and weaning practices of ancient human populations 
can provide information on the health, cultural traits, and reproduction of these populations. 
Breast milk provides various antibodies as well as nutrition to infants, and is important for 
subadult survival [IJ[2] . Breastfeeding practices are closely related to the growth of subadults and 
overall health of a population [3-5J. The type of subsistence activities, social constructs, diet and 
various cultural factors affect breastfeeding practices [6"8], and the length of the breastfeeding 



period is one of the most important determinants of the fertility of a population J9J[T0]. Shorter 
breastfeeding periods tend to result in shorter birth intervals, and, in turn, higher fertility because 
breastfeeding can delay the resumption of ovulation |llfll4) . Furthermore, it is supposed that 
humans are weaned earlier than the other great apes, and understanding evolutionary changes 
in weaning practices is of great interest |15H2Q| . 

Stable isotope analysis of bone collagen is useful for reconstructing the dietary habits of 
ancient people, and it has also been used to reconstruct breastfeeding and weaning practices 
of archaeological populations [21H25J . Nitrogen isotope ratios (£ 15 N values) of body proteins 
primarily reflect dietary protein isotope ratios [26, 27|. Prior to and immediately after birth, 
<5 15 N values of infants are the same as those of their mothers [28J. After birth, infants who are 
exclusively breastfed show 2-3%o higher <5 15 N values than their mothers |22|,I29] because of the 
trophic level effect |30H32| . Subadult <5 15 N values decrease after the introduction of supplemen- 
tary foods, and gradually approach the values found in adult bone collagen. It is possible to 
reconstruct infant feeding practices of an archaeological populations by combining <5 15 N values 
and physically estimated ages at death of subadults of different ages [5|I33|. 

However, in previous isotopic studies, weaning ages have been subjectively estimated from 
visual assessments of detectable changes in subadult bone collagen (5 15 N values. To overcome 
these difficulties, attempts have been made to simulate changes in (5 15 N values of subadult bone 
collagen in two pioneering studies. Schurr [34] used exponential functions to describe changes in 
<5 15 N values and estimate the age at the start of weaning. Millard [35j suggested that the model 
proposed by Schurr |34J suffered from a number of difficulties, and proposed an alternative model 
that further included a nitrogen mass balance and the age at the end of weaning. However, both 
models still suffer from the following three problems. 

1. The subadult bone collagen turnover rates are not fully considered. The bone collagen 
turnover rate is high in early infancy |36[I37|. but it decreases over the course of subadult 
growth [38, 39j . If not corrected, the lower bone collagen turnover rates at higher ages 
would generate significant discrepancies between the visible changes in bone <5 15 N values 
and actual weaning ages. 

2. Some parameters used to describe changes in <5 15 N values are determined arbitrarily. Two 
parameters, 15 N-enrichment from maternal to infant tissues and the <5 15 N values in weaning 
foods, could vary among different individuals and populations; therfore, they should be 
considered as variables in addition to the weaning ages. First, it has been reported that 15 N- 
enrichment varies to some extent in modern infant-mother pairs (between 1.7%o and 2.8%o, 
n = 7: |29j) and in archaeological populations (between 0.5%o and 4.4%o, n = 25: [40J). 
Second, it is possible that <5 15 N values of materials used in weaning foods were different 
than those used in adult foods 14111421. 



3. The results are represented as point estimates without either probabilities or confidence 
intervals. The probabilities of the weaning parameters should be calculated to evaluate the 
validity of the computation results. 



The objective of this study is to develop a model for analyzing cross-sectional <5 15 N data of 
subadult bone collagen in archaeological skeletal populations. The model is programmed in R 
language, which is a free software environment for statistical computing and graphics [43j. The 
model has the following three important features that are not present in the previous models: 

1. The subadult bone collagen turnover rate is estimated anew and incorporated in the equa- 
tions. 

2. The enrichment factor and <5 15 N values of weaning foods are included as target parameters 
to be estimated. 

3. Using a framework of approximate Bayesian computation (ABC) allows researchers to 
calculate the probabilities and credible intervals of the weaning parameters. 



Subadult bone collagen turnover rate 

Temporal changes in the bone collagen turnover rate must be considered to estimate a precise 
weaning ages from an observed isotope ratio. Bone collagen is laid down during childhood because 
of bone modeling, which is a formative process primarily associated with skeletal growth, and is 
replaced throughout life by bone remodeling, which is a coupled resorptive and formative process 
that does not change the quantity of bone [44, 45j. As indicated in Figure [H turnover refers to 
the proportion of newly synthesized bone collagen to the total bone collagen during modeling 
and remodeling over a unit of time. When the turnover rate is high enough (i.e., >1.0 per unit 
time), bone collagen at a specific age consist only of newly synthesized collagen, and the isotope 
ratio will immediately change with dietary changes. When the turnover rate is lower (i.e., <1.0), 
the bone collagen consists not only of newly synthesized but also previously synthesized collagen, 
the isotope ratio reflects recent and past dietary intakes. 
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Figure 1. Schematic illustration of the bone turnover process. The <5 15 N value for 
bone collagen at the unit time age of t years is represented as 5 15 A^b one (i). 



Although temporal changes in turnover rates of subadult bone minerals and collagen have 
been estimated by analyzing the uptake of 90 Sr fallout [37, 46J and bomb- 14 C [38J, respectively, 
the estimates produced from these bulk cross-sectional studies were not necessarily precise. Some 
of the estimates were not based on direct measurements in subadults but on extrapolations from 
results for adults. In addition, the assumptions made about the dietary intake of tracers in these 
subadults were simplistic and ignore individual variation. In the present study, we calculated 
turnover rates from bone metabolism mechanisms at the tissue level so that more precise subadult 
bone turnover rates could be estimated. 

Turnover rates of mineral and organic phases should differ because the mineralization process 
is much slower than the synthesis of the organic matrix. Bone is a composite material, and is 
mainly made of a calcified organic matrix [45 f 47j. The microstructure of bone material consists of 
assembled collagen fibrils forming the organic phase, and tiny mineral particles reinforcing them 
|44j . Two coupled processes are responsible for bone remodeling. The resorptive process involves 
osteoclasts dissolving the mineral phase by creating a low pH environment around the bone 
surface, and then producing a lysosomal protease to degrade the organic matrix |48| . The next 
formative process inolves osteoblasts replacing the organic matrix and rapidly mineralizing it to 
up to 70% of full mineralization capacity within a few days (primary mineralization) , the residual 
30% of the mineralization occurring gradually over several years (secondary mineralization) |44[ 
149] . The mineralization process has been formulated as a mineralization law |49| . Bone modeling 
occurs with a similar formative process as bone remodeling but with the resorption of the bone 
cartilage template instead of mineralized old bone [33] . Since the growth |50] and replacement |51| 
(i.e., turnover, consisting of the modeling and remodeling processes) of bone minerals at the tissue 
level have been well documented, the turnover of bone collagen can be estimated by correcting 
the mineralization delay 



Approximate Bayesian computation 

ABC is a modern approach in Bayesian inference that allows posterior distributions to be evalu- 
ated when it is difficult to calculate the likelihood function, which describes probabilities under 
given parameters. Various ABC methods have been applied in diverse fields such as population 
genetics, evolutionary biology, ecology, and epidemiology [52H54] . 

A general ABC algorithm takes a given observation x and repeat the following three steps 
until J points have been accepted: 

1. Draw the candidate parameter 9j from the prior distribution n(0). 

2. Simulate dataset Xj using 9j and the model. 

3. Accept 0j if p(x,Xj) < a, and otherwise reject Qj. 

Here p(-) is a function measuring the distance between simulated and observed data points, a is a 
fixed tolerance for the "closeness" of simulated and observed data, and x, Xj, and 9 may be vector 



values. If p(-) measures appropriate distances and tolerance is sufficiently small, the accepted 
parameters reasonably approximate the posterior distributions. This is a rejection sampling 
algorithm, which is the simplest ABC procedure. 

Although ABC has proved to be a flexible and powerful approach for evaluating posterior 
distributions, its major drawback is its inefficiency. Acceptance rates in the simple rejection sam- 
pling described above can be very low, especially when the posterior is a long way from the prior, 
which wastes computing time. Several algorithms have been proposed to increase the sampling 
efficiency, by introducing weighting with regression analysis [55,56], Markov chain Monte Carlo 
sampling [57], and sequential Monte Carlo (SMC) sampling |58^60| . We used SMC sampling 
with corrected partial rejection control proposed by Sisson et al. |59] because this method could 
be implemented more quickly and simply in our model in the R software environment. SMC 
sampling is characterized by successively decreasing the tolerance, and weighted resampling from 
the previous parameter population. 

Materials and Methods 

Estimating subadult bone collagen turnover rates 

In this study, bone collagen turnover rates in subadults were calculated from the modeling |50| 
and remodeling |51] rates for cancellous bone minerals, and the mineralization law for the bone 
organic matrix [49] . "Turnover" is defined as the aggregated effects of bone modeling (i.e., the 
addition of bone tissue by skeletal growth) and remodeling (i.e., the replacement of existing bone 
tissues). First, following Leggett et al. [51] . the bone mineral turnover rate T m j n [t] over one unit of 
time (i.e., one year from t—1 to t, Equation 4) in childhood was calculated using the functions that 
describe the temporal change in bone mineral mass C(t) |50| (Equation 1) and the remodeling rate 
j(t) [51] (Equation 2). Next, the bone collagen turnover rates T co /[i] over one unit of time from 
t — 1 to t years (Equation 5) were calculated sequentially with T m i n [t] and the mineralization law, 
X(i), which described the bone collagen mineralization process [49] . The mineralization law was 
derived from Ruffoni et al. [49] . and represents the rate of mineralization of the collagen portion 
at the ith year after the collagen matrix was formed (Equation 3). Finally, the resulting discrete 
turnover rates were coerced into a quartic polynomial (QP) formula (Equation 6). Turnover 
rates at ages less than one year were extrapolated from the QP function. 

Basic functions to describe the temporal changes of bone mineral were derived from several 
previous studies. Following Mitchell et al. [50j, the bone mineral mass C(t) at age of t years was 
represented as: 

C(t) = 28.0 + 86.828t - 16.5105£ 2 + 1.5625t 3 - 0.04114t 4 (0 > t > 20) (Equation 1). 
This equation represents the modeling process of bone turnover. On the other hand, the remod- 
eling rate j(t) at age of t years was represented as follows: 

1. j(f) = ^&f , when t < 1.5, and 



2. j(t) = 0.975e~ am , when t > 1.5 (Equation 2). 



This equation was obtained from Leggett et al. [51J and was derived from direct histological 

observations of subadult rib bone formation and resorption performed by Frost |61j . Note that 

a term for the radioactive decay of 90 Sr (0.025 per year), included in the original equations, was 

excluded from our equations. Following Ruffoni et al. |49| . the mineralization law X(i), which 

describes the rate of mineralized collagen portion at ith. years after the collagen matrix was 

formed, was set as: 

1+4 1+4 

X(i) = c\ _ L n + C2 _/ 2 (Equation 3). 

»1 i 2 

In this study, c\, i\, c<i and i<i were set as 18/23, 1/300, 25/92, and 5, respectively. Equation 
3 corresponds to over 70% primary mineralization in a few days and protracted secondary min- 
eralization of up to 100% in about 20 years, values that have been given in several previous 

studies @umgE2iiS3]. 

Temporal changes in the bone mineral turnover rate can be represented using these functions. 
Put simply, the bone mineral turnover rate, T m i n [t], over one unit of time from t — 1 to t years, 
was represented as follows: 

T m in[t] = - gft) + ft-1 r y( x ) dx cW (Equation 4). 

The former and latter terms in the function indicate the effects of bone modeling and remodeling, 
respectively. Using the bone collagen turnover rate, T co ;[t], over one unit of time from t — 1 to t 
years, T m i n [t] can also be represented as follows: 

1- T min [t] = Tcoi [t]AX[l], when t = 1, and 

2. T min [t] = T col [t]AX[l] + E%\(T C oi[j]^AX[t + 1 - j]), when t > 2 (Equation 5). 

The former and latter terms in the second function shown in Equation 5 indicate the effects of 
turnover delay in the bone mineral for the intended unit of time (i.e., t — 1 to t years) and the 
aggregated effects of the delay for the former unit times (i.e., to 1 year, 1 to 2 years, .., and 
t — 2 to t — 1 years). The turnover rate over one unit of time (i.e., one year) from t — 1 to t 
years can be sequentially calculated using Equation 4 and 5. The resulting discrete turnover 
rates were coerced to a quartic polynomial formula using the nls function in R, in accordance 
with the quartic formula for bone mineral mass (i.e., Equation 1). The formula is represented as 
follows: 
T col [t] = 1.778 - 0.412U + 0.05029t 2 - 0.002756£ 3 + 0.0005325t 4 (Equation 6). 

Changes in 5 15 N values of diet and bone collagen 

Following Millard [35], the 5 l5 N value of newly synthesized collagen at a given age of t years was 
defined by four parameters, the ages at the start (t\) and end (£2) of weaning, enrichment factor 
between the infant and mother (E), and 5 15 N value of collagen synthesized entirely from weaning 
foods (5 15 N wn f ooc [) (Equation 7 and 8). The <5 15 N value of newly synthesized collagen equals the 
sum of the 5 15 N value of the mothers tissue and enrichment factor before weaning (t < t\), 
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which changes exponentially during weaning (t\ < t < £2)) an d equals the collagen <5 15 N value 
that fully reflects the consumption of supplementary food (t > £2)- Then, the incorporation 
of newly synthesized collagen and replacement of existing collagen in bone are simulated in 
over each successive unit time using the estimated turnover rate for bones (Equation 6 and 9). 
As most isotopic studies on weaning have focused on rib bones, because of their assumed fast 
turnover [64] and relatively trivial importance in morphological studies, the rate incorporated 
into the present model was that of cancellous bones. Although the rib bones that were sampled 
would have contained cortical parts, the relatively high surface to volume ratio in ribs would 
have resulted in a high proportion of cancellous parts and only thin cortical parts, making the 
turnover rate comparable to that of cancellous bones |64] . Although one unit of time consists of 
one year, adjustments from the last unit of time enables simulated <5 15 N values to be calculated 
for each individual in the dataset (Equation 10 and 11). Simulated <5 15 N values, 5 15 N( )one , for 
each individual can be calculated under the given weaning parameters (t\, t2, E and 5 15 N wn f ooc i) 
using the model described above. The most appropriate weaning parameters can be estimated by 
minimizing the mean least square distance between the observed and resultant simulated change 
in bone collagen <5 15 N values. 

Following Millard [35], the <5 15 N values for newly synthesized collagen <5 15 N new (t) at a given 
age of t years are given by the following equation: 

5 15 N new {t) = {l-p(t)){6 15 N mother + E) +p(t)5 15 N wnfood (Equation 7). 

The proportion of non-milk protein in the total dietary protein intake at age of t is represented as 
p(t). The <5 15 N value for the mother's milk is described as <5 15 iV moi /j er . + E, using the 5 15 N value 
for the mothers tissue, S 15 N mot h er (approximated by the mean 5 15 N value for adult females), 
and a 15 N enrichment factor for the transfer from the maternal to infant tissue, E. The <5 15 N 
value for collagen synthesized from non-milk foods is represented as 5 15 N wn f 00 d. We considered 
5 15 N wn f ooc i to be a variable because children in the past could have eaten supplementary foods 
with different <5 15 N values from the adult mean <5 15 N values. This value has been approximated 
in previous studies as the mean <5 15 N value for the adults. 

The proportion of non-milk protein in the total dietary protein intake is assumed, in our 
model, to increase exponentially. The relative proportion of non-milk protein at the age of t 
years, pit), is described as follows: 

1. pit) = 0, when t < t\ (breast milk only), 

2. p{t) = { t ~_l ) 2 , when t\ < t < £2 (during the weaning process), and 

3. pit) = 1, when t > t<i (no breast milk), (Equation 8), 

where the ages at the start and end of weaning are represented as t\ and t2, respectively. Equation 
8 was derived from a model proposed by Millard |35| and represents slow initial weaning and 
rapid final weaning. In the original model, four forms (linear, parabolic, reverse parabolic, 
and sigmoid) of dietary change were applied to condition 2 in Equation 8. Although the form 
of dietary change during weaning can be selected in the WARN package, we used only the 



parabolic form because Millard |35j used a parabolic weaning pattern to model the changes 
in <5 15 N in archaeological datasets. This seems to be a reasonable assumption as not only the 
amount of milk protein consumed decreases during the weaning process but also the proportion 
of milk protein consumed also decreases because of the increasing total dietary intake in growing 
subadults. 

The 8 15 N value for bone collagen at the age of t years, 5 15 N bone (t), is calculated as follows: 
S 15 N bone (t) = S 15 N bone (t - 1)(1 - T col [t}) + //_! 5 15 N new (t){x)dxT co i[t] (Equation 9). 
The former and latter parts of the equation represent the remaining and the newly synthesized 
portion, respectively, of the bone collagen over one unit of time from t — 1 to t years. Extending 
equation 9, the <5 15 N value for bone collagen at the age of t + a, i.e., a being a part of one year 
from the unit time point t (0 < a < 1), is represented as: 

S 15 N bone {t + a)= 5 15 N bone (t){l - T col [t + a}) + // +a 5 15 N new (t){x)dxT col [t + a] (Equation 10). 
In equation 10, T co /[T] is the bone collagen turnover rate over a year from t to t + a, given by: 

T col [t + a} = T col [t + l] jiZl2lw2 ( E q uation 11)- 

The bone collagen <5 15 N values for each unit of time (one year) can be calculated sequentially, as 
reference values, using Equation 9 under the given parameters. The <5 15 N values that correspond 
to the observed ages for the samples can then be calculated from the reference values and Equa- 
tion 10. The initial bone collagen values at year of age, 5 15 N bone (0), were approximated using 
the mean e) 15 N value of adult females, because the 5 15 N value for infant tissue is assumed to be 
the same as to that of the mother |29| . Theoretical <5 15 N values for the age of each individual in 
the observed dataset can be calculated using Equation 10. 

In our model, the differences between the individuals are evaluated by calculating mean square 
distance, D, between the observed and simulated <5 15 N values. Put simply, point estimates of 
the parameters with minimized D can be calculated by solving the optimization problem (the 
application of the optimization problem to palaeo dietary reconstructions has been described 
by Little and Little |65j). These represent point estimates under the framework of maximum 
likelihood estimates (MLE) . Although the point estimates do not provide information on the error 
ranges, they will be used later in the SMC sampling procedures; therefore, optimized values for 
weaning parameters under the MLE framework were calculated. We used the optim function 
in R to obtain the optimized parameter value, 9 op t, and its resultant minimum mean square 
distance, D opt . 

Incorporation of ABC 

To obtain posterior probabilities of the estimated parameters, fitting calculations between the 
observed and simulated data are performed under the ABC framework with SMC sampling pro- 
posed by Sisson et al. [59j. Using the ABC framework, a number of weaning parameter sets that 
give well- fitted <5 15 N values were sampled and assumed to represent the posterior distributions 
of the parameters. After applying the ABC procedure, posterior distributions were smoothed 
using the kernel density estimation |66j . and joint probabilities for weaning ages (t± and £2) and 



marginal probabilities for E and 5 15 N wn f ooc i were calculated. In the density estimation, poste- 
rior probabilities were calculated to one decimal places for discrete parameter categories because 
strictly implementing the density estimation as a continuous distribution requires advanced nu- 
merical analysis techniques. 

SMC sampling is characterized by a successive reduction in tolerance and a weighted re- 
sampling from the previous parameter population, called a "particle". Particles of preliminary 
simulations are used to calculate the next set of parameter vectors, to generate simulated data 
within a certain distance D from the observed data. The particles are then repeatedly resam- 
pled (according to a weighting scheme that considers the prior distributions), perturbed (using 
a transition kernel), and judged (on the basis of a successively decreasing tolerance). The par- 
ticles after this iterative process finally approximate a sample of the posterior distribution of 
the parameters. In particular, the partial rejection control procedures prune away parameters 
that have minimal impacts on the final estimation in the parameter weighting step in the earlier 
stages of the tolerance reduction, and this increases the sampling efficiency |67j . 

To adopt the ABC framework, we added individual error terms t% in Equation 10 as follows: 
S l5 N bone (t + a) = S 15 N bone (t)(l - T col [t + a}) + f f t+a 8 15 N new (x)dxT col [t + a] + e t (Equation 12). 
These errors were independently sampled from the normal distribution with mean of 0.0 and SD 
of a, and individually assigned to simulated <5 15 N values. By considering this individual error 
term, parameters that result in D values smaller than D op t can be generated, which represent 
more plausible estimates for the measured data. In the ABC framework, D values are calculated 
using randomly generated parameters from the prior distributions, then the parameters that 
result in D values smaller than D op t become the posterior distributions. 

The sequential Monte Carlo algorithm in our model proceeds as follows (see Sisson et al. |59| 
for more details): 

1. Set prior distributions tt(-) for the parameters and the number of particles j in one pop- 
ulation. Calculate the final tolerance ax (= Dopt) under the MLE framework and set 
decreasing tolerances. Set the population indicator k — 1 (initialization). 

2. Set the particle indicator j — 1 (initialization). 

(a) If k = 1, independently sample 9** from the prior distribution n(0). If k > 1, sample 

(i) (i) 

9* from the previous population 6\_^ with weights Wu_ 1: and perturb the particle to 

9** with transition kernel 0. Simulate the change in the <5 15 N value S 15 N** w (t) with 

9** using equation 12. If D** — D op t > cufc, 9** are rejected and then repeat procedure 

2(a). 

(b) Set the indicators as follows: 



(if fc>l). 



mJ) /]** 




W { k j) = 1 (if k = 


- 1), and 


W^ - 
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If j < J, increment j = j + 1 and go to procedure 2(a). 

3. Normalize the weights so that: 

E/=i wjp = i- 

If the requirements for an effective sample size ESS are not met such as: 

ESS = — t < o, 

E^>F 2 ' 

(7) (7) 

sample with replacement, the particles 9\ , with weights W^ to obtain a new population 

e[ J \ and set weights W^ j) = j. 

4. If A; < i^, increment k = k + 1 and go to procedure 2. 

Default prior distributions 7r(-) were set as normal distributions with default means of {0.5, 3.0, 
1.9, S 15 N mother , and 0.0} and SDs of {3.0, 3.0, 0.9, 3.0, and 1.0} for t u t 2 , E, S 15 N wnfood , 
and a, respectively. The mean weaning age was obtained from values recommended by modern 
pediatricians and the biologically expected ages |68j . The mean and standard deviation of the 
enrichment factor E was obtained from the values reported by Waters- Rist and Katzenberg |40j . 
The hyper parameter for the individual error term a was used as an absolute value in the 
calculation. The default number of particles J was 10000. Decreasing tolerances a^ were set as 
D op t + {2, 1, 0.5, 0.25, 0.125, 0.0625, 0} and, therfore, the number of populations K — 7. The 
transition kernel (j) was set to be a normal distribution with a mean of 0.0 and SD of 0.1. 

Results and Discussion 

Subadult bone collagen turnover rate 

The calculated turnover rates are shown in Table Q] and Figure [2J The turnover rate of bone 
collagen was estimated to be larger than that of bone mineral until an individual reaches their 
late teens, and to decrease over the course of subadult growth. 

The integrated bone collagen turnover rate from 0.0 to 1.0 years of age was estimated to 
be 1.588, and the estimated bone collagen turnover rate was higher than 1.000 per year by two 
years of age (see Table [1]). The integrated turnover rate from 0.0 years of age reached 0.966 at 
0.60 years of age, suggesting that it takes 31 weeks for infants to fully reflect post-birth dietary 
<5 15 N signals. Tracer intake and biochemical marker studies have shown that the bone mineral 
and collagen turnover rates are high in the first few years of life (i.e., >1.0 per year) |36 [ l37 1 l39] . 
which is consistent with our results (see Tableland Figure [2]). However, temporal changes in 
the bone collagen turnover rate after infancy and before adulthood have never been estimated 
directly and continuously, and the present study allowed them to be estimated. An isotopic 
study on an archaeological infant of a known age has suggested that infant rib bone collagen can 
fully reflect post-birth dietary 15 N input, in an extreme case, in only five to six weeks [69], but 
this is estimated to take 31 weeks from our results. Our study allows typical temporal changes 
to be estimated, but the bone collagen turnover rate in subadulthood probably varies. 
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Table 1. Estimated temporal changes in turnover rates for bone minerals and 
collagen. 



Age 




Turnover rate 




From 


To 


Mineral 


Collagen 


Collagen (QP) 





1 


1.217 


1.474 


1.413 


1 


2 


0.908 


1.059 


1.134 


2 


3 


0.786 


0.892 


0.924 


3 


4 


0.700 


0.776 


0.771 


4 


5 


0.629 


0.682 


0.664 


5 


6 


0.571 


0.611 


0.590 


6 


7 


0.527 


0.558 


0.540 


7 


8 


0.492 


0.520 


0.507 


8 


9 


0.462 


0.489 


0.483 


9 


10 


0.434 


0.461 


0.463 


10 


11 


0.407 


0.432 


0.441 


11 


12 


0.378 


0.402 


0.416 


12 


13 


0.349 


0.370 


0.386 


13 


14 


0.319 


0.337 


0.349 


14 


15 


0.289 


0.302 


0.306 


15 


16 


0.258 


0.267 


0.260 


16 


17 


0.227 


0.231 


0.213 


17 


18 


0.194 


0.193 


0.171 


18 


19 


0.158 


0.151 


0.139 


19 


20 


0.118 


0.104 


0.124 



QP: calculated from the QP function. 



The integrated bone collagen turnover rate from 19.0 to 20.0 years of age was estimated to 
be 0.130 per year in our study (see Figure [2]), which is a little higher than that proposed by 
Stenhouse and Baxter (10.4 ± 2.7% during adulthood, [70]) and Hedges et al. (9.7% and 4.1% 
for 20- year-old male and female femora, respectively, [38J ) . Although the type of bone sampled 
by Stenhouse and Baxter |70j is not stated, differences between the turnover rates in different 
bone types could cause these different results. The turnover rates are higher in bones with greater 
surface to volume ratios than those in bones with smaller ratios |64J. Ribs, which were target 
bones in our study, have relatively high proportions of cancellous and thin cortical parts, whereas 
femur analyzed by Hedges et al. |38| has a lower proportion of cancellous and thick cortical parts, 
although there are slight differences, the overall trend of the temporal changes in bone turnover 
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Figure 2. Estimated temporal changes in bone mineral and collagen turnover 
rates. Turnover rates of bone minerals and collagen are represented as discrete values, and 
that of collagen is fitted to QP plotted against age. 

rates in this study is consistent with previous estimates. 

The implemented model 

The model developed in the present study is implemented as the R package WARN (Weaning Age 
Reconstruction with Nitrogen isotope analysis). Credible intervals can be calculated for a given 
parameter range using the WARN package. Images of the results calculated using the package 
are shown in Figure [3l Application of this model to previously reported skeletal populations and 
meta-analysis of the results will be reported elsewhere in the near future. 

Although it is desirable to test the model validity, the absence of proper test data means this 
is not possible. Archaeological skeletal populations cannot be tested because the true weaning 
ages are usually unknowable, and historical literature, if any, describing breastfeeding practices 
at the time period when the population lived sometimes differs from actual practices (e.g. [69U72J ; 
see also |41 [ l73 [ l74]). Since the model presented here was intended for human subadult bones, 
conducting an experimental study was difficult, and hair, nail, and other tissues were not suitable 
for analysis because they have different turnover rates than bone collagen. Experimental studies 
of animals would not be appropriate because human growth patterns are unique among mammals 
|18p75| ; therefore the nitrogen mass balance in human subadults would probably be different from 
that in other animals. 

There are two caveats to consider before applying the model presented here. First, the 
present model is intended for bones with relatively high turnover rates, such as cancellous bones 
or ribs. Although WARN can be applied equally to isotopic data from bones with relatively 
low surface to volume ratios (e.g., limbs, cranium, and mandible), attention to this aspect is 
required for more precise analysis. Second, the WARN approach will always attempt to fit a 
model, even if the subadult <5 15 N values do not indicate breastfeeding and weaning signals. If 
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Figure 3. An example of the results of applying WARN model using the 
Spitalfields population as a case study. (A) Modeled temporal changes in the <5 15 N values 
by subadult age calculated from the reconstructed MDEs. Mean and SD ranges for adult 
females and all adults are indicated with open circles and crosses, respectively. (B) Contour 
lines show the posterior probability for the combination of weaning ages. The target ranges for 
t\ and ti are 0.0-1.2 years and 1.2-2.0 years of age, respectively, and the calculated joint 
probability for the ranges is 0.956. (C) Distribution of posterior probabilities for the 
15 N-enrichment from maternal to infant tissues. The target range is 1.6— 2.4%o, and the 
calculated marginal probability for the range is 0.967. (D) Distribution of posterior 
probabilities for the <5 15 N values for collagen synthesized entirely from weaning foods. The 
target range is 12.4— 13.0%o, and the calculated marginal probability for the range is 0.975. 
Subadult ages and bone collagen 5 15 N values were obtained from Nitsch et al. |69|I71|. 



researchers cannot find patterns of isotopic changes by visually inspecting the data, they are 
urged to examine their data carefully before applying the model, for example, for a biased age 
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distribution or high isotopic variability in subadults. Although the estimated turnover rate and 
model developed can be further improved, in this study, we propose a framework for objectively 
and quantitatively analyzing and interpreting subadult bone collagen e) 15 N values. A precise 
reconstruction of past breastfeeding and weaning practices over a wide range of time periods and 
geographic regions could make it possible to understand this unique feature of human life history 
and cultural diversity in infant feeding practices |15H20j . 
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